*******************************************************************************
* .do-file for section on housing robustness
*******************************************************************************

		clear all
		set maxvar 6000
		cd "$scfwkddata"
		
		use dataset_analysis, clear
		
*************************************************************
************* prepare for multiple impuation ****************
*************************************************************		

		tostring Y1, replace
		g imp = substr(Y1,-1,.)
		
		tostring wave, replace
		g Y1_uni = wave+Y1
		
		destring Y1 imp wave Y1_uni, replace
		
		
*************************************************
********** final variable adjustment ************
*************************************************
		
		// dummy indicator for home ownership
		gen homeowner = .
		* lives in house/apartment and owns all, part, or lives in townhouse association
		replace homeowner = 1 if hometype == 3 & (owntype_oth == 1 | owntype_oth == 8 | owntype_oth == 5)  
		* lives in mobile home and own both site and home, or either of the two 
		replace homeowner = 1 if hometype == 2 & (owntype_mob == 1 | owntype_oth == 2 | owntype_oth == 3)
		* lives on farm or ranch, and owns all or owns part
		replace homeowner = 1 if hometype == 1 & (owntype_farm == 1 | owntype_farm == 2)
		* all others
		replace homeowner = 0 if homeowner == .
		
		// further housing variables: housing-wealth-to-income-ratio
		gen HI = housing / (totinc_nocapoth ) if (totinc_nocapoth ) > 0
		sum HI, det
		replace HI = . if HI > r(p99) 
		
		// only compare singles
		drop if single == .
		
		// set as survey data
		svyset Y1 [pweight = wgt_rv] 
		
		**********************************************************************************
		****** First exercise: LC Profile of Housing Patterns (Singles) ******************
		**********************************************************************************

		// homeownership rate single men vs single women  
		gen ci_low  = . 
		gen ci_high = . 
		gen meanhr    = . 
		
		foreach x of numlist 0/1 {
		foreach num of numlist 30(1)65 {
		quietly svy: mean homeowner if ageav == `num' & fam_type == `x'
		mat A = r(table) 
		replace  meanhr = A[rownumb(A, "b"),1 ] if ageav == `num' & fam_type == `x'
		replace ci_low  = A[rownumb(A, "ll"),1 ] if ageav == `num' & fam_type == `x'
		replace ci_high = A[rownumb(A, "ul"),1 ] if ageav == `num' & fam_type == `x'
		}
		}

		
		sort ageav		
		tw (line ci_low ageav  if fam_type == 0, lp(solid) lc(black*0.15) lw(medium)) ///
		(line meanhr  ageav    if fam_type == 0, lp(solid) lc(black) lw(medium)) ///
		(line ci_high ageav    if fam_type == 0, lp(solid) lc(black*0.15) lw(medium)) ///
		(line ci_high ageav    if fam_type == 1, lp(solid) lc(black*0.15) lw(medium)) ///
		(line ci_low  ageav    if fam_type == 1, lp(solid) lc(black*0.15) lw(medium)) ///
		(line  meanhr ageav    if fam_type == 1, lp(dash)  lc(black) lw(medium)),  ///
		xtitle("Age", siz(vlarge)) ytitle("Homeownership Rate", size(vlarge))  ///
		legend(order(2 "Men" 6 "Women") size(vlarge) ///
		position(6) col(4) region(lcolor(white)) stack) ///
		graphregion(color(white)) xlabel(30(10)65,labsize(vlarge)) ylabel(0.1(.2)0.7,labsize(vlarge))
		graph export "$resultpath/singles_homeownershiprate.eps", replace
		graph export "$resultpath/singles_homeownershiprate.pdf", replace
		
		
		// gross housing wealth single men vs single women
		replace ci_low  = . 
		replace ci_high = . 
		gen meanhous    = . 
		
		foreach x of numlist 0/1 {
		foreach num of numlist 30(1)65 {
		quietly svy: mean housing if ageav == `num' & fam_type == `x'
		mat A = r(table) 
		replace meanhous = A[rownumb(A, "b"),1 ] if ageav == `num' & fam_type == `x'
		replace ci_low   = A[rownumb(A, "ll"),1 ] if ageav == `num' & fam_type == `x'
		replace ci_high  = A[rownumb(A, "ul"),1 ] if ageav == `num' & fam_type == `x'
		}
		}
		replace meanhous = meanhous/1000
		replace ci_low = ci_low/1000
		replace ci_high = ci_high/1000
		
		sort ageav
		
		tw (line ci_low ageav  if fam_type == 0, lp(solid) lc(black*0.15) lw(medium)) ///
		(line meanhous ageav   if fam_type == 0, lp(solid) lc(black) lw(medium)) ///
		(line ci_high  ageav   if fam_type == 0, lp(solid) lc(black*0.15) lw(medium)) ///
		(line ci_high  ageav   if fam_type == 1, lp(solid) lc(black*0.15) lw(medium)) ///
		(line ci_low   ageav   if fam_type == 1, lp(solid) lc(black*0.15) lw(medium)) ///
		(line meanhous ageav   if fam_type == 1, lp(dash)  lc(black) lw(medium)),  ///
		xtitle("Age", siz(vlarge)) ytitle("Hous. Wealth in (000s) 2007\$", size(large))  ///
		legend(order(2 "Men" 6 "Women") size(vlarge) ///
		position(6) col(4) region(lcolor(white)) stack) ///
		graphregion(color(white)) xlabel(30(10)65,labsize(vlarge)) ylabel(0(50)200,labsize(vlarge))
		graph export "$resultpath/singles_housingwealth.eps", replace
		graph export "$resultpath/singles_housingwealth.pdf", replace

		
		// Housing-to-income ratio single men vs single women
		replace ci_low  = . 
		replace ci_high = . 
		gen meanHI      = . 
		
		foreach x of numlist 0/1 {
		foreach num of numlist 30(1)65 {
		quietly svy: mean HI if ageav == `num' & fam_type == `x'
		mat A = r(table) 
		replace meanHI   = A[rownumb(A, "b"),1 ] if ageav == `num' & fam_type == `x'
		replace ci_low   = A[rownumb(A, "ll"),1 ] if ageav == `num' & fam_type == `x'
		replace ci_high  = A[rownumb(A, "ul"),1 ] if ageav == `num' & fam_type == `x'
		}
		}
		
		sort ageav
		tw (line ci_low ageav  if fam_type == 0, lp(solid) lc(black*0.15) lw(medium)) ///
		(line meanHI ageav   if fam_type == 0, lp(solid) lc(black) lw(medium)) ///
		(line ci_high  ageav   if fam_type == 0, lp(solid) lc(black*0.15) lw(medium)) ///
		(line ci_high  ageav   if fam_type == 1, lp(solid) lc(black*0.15) lw(medium)) ///
		(line ci_low   ageav   if fam_type == 1, lp(solid) lc(black*0.15) lw(medium)) ///
		(line meanHI ageav   if fam_type == 1, lp(dash) lc(black) lw(medium)),  ///
		xtitle("Age", siz(vlarge)) ytitle("Hous. Wealth-to-Income Ratio", size(large))  ///
		legend(order(2 "Men" 6 "Women") size(vlarge) ///
		position(6) col(4) region(lcolor(white)) stack) ///
		graphregion(color(white)) xlabel(30(10)65,labsize(vlarge)) ylabel(,labsize(vlarge))
		graph export "$resultpath/singles_housingwi.eps", replace
		graph export "$resultpath/singles_housingwi.pdf", replace
			
		
		*****************************************************************************
		****** Second Exercise: Raw data of Equity Share  by Housing Tenure *********
		*****************************************************************************
		
		replace ci_low   = . 
		replace ci_high  = . 
		gen mshare    = . 
		foreach x of numlist 0/1 {
			foreach y of numlist 0/1 {
		foreach num of numlist 30(1)65 {
		quietly svy: mean share if ageav == `num' & single  == `x' & homeowner == `y'
		mat A = r(table) 
		replace mshare   = A[rownumb(A, "b"),1 ] if ageav == `num' & single  == `x' & homeowner == `y'
		replace ci_low  = A[rownumb(A, "ll"),1 ] if ageav == `num' & single  == `x' & homeowner == `y'
		replace ci_high = A[rownumb(A, "ul"),1 ] if ageav == `num' & single  == `x' & homeowner == `y'
		}
		}
		}
		
		gen std_dev = (ci_high-mshare) / 1.96
		
		preserve
		keep mshare ageav single homeowner ci_high ci_low
		duplicates drop
		sort ageav
		tw (line ci_high ageav if homeowner == 0 & single == 1,          lc(black*0.15) lw(medium)) ///
		   (line ci_low  ageav if homeowner == 0 & single == 1,          lc(black*0.15) lw(medium)) ///
		   (line mshare  ageav if homeowner == 0 & single == 1, lp(dot)  lc(black) 		lw(medium)) ///
		   (line ci_high ageav if homeowner == 1 & single == 1, 		 lc(black*0.15) lw(medium)) ///
		   (line ci_low  ageav if homeowner == 1 & single == 1, 		 lc(black*0.15) lw(medium)) ///
		   (line mshare  ageav if homeowner == 1 & single == 1, lp(dash) lc(black) 		lw(medium)) ///
		   (line ci_high ageav if homeowner == 0 & single == 0, 		 lc(black*0.15) lw(medium)) ///
		   (line ci_low  ageav if homeowner == 0 & single == 0, 		 lc(black*0.15) lw(medium)) ///
		   (line mshare  ageav if homeowner == 0 & single == 0, 		 lc(black) 		lw(medium)) ///
		   (line ci_high ageav if homeowner == 1 & single == 0, 		 lc(black*0.15) lw(medium)) ///
		   (line ci_low  ageav if homeowner == 1 & single == 0, 	     lc(black*0.15) lw(medium)) ///
		   (line mshare  ageav if homeowner == 1 & single == 0, lp(dash_dot) lc(black)  lw(medium)), ///
		xtitle("Age", siz(vlarge)) ytitle("Risky Share in %", size(vlarge))  ///
		legend( order(9  "Men (renter)" 12 "Men (owner)" 3 "Women (renter)" 6 "Women (owner)" )size(vlarge) position(6) /// 
		col(2) region(lcolor(white)) stack) ///
		graphregion(color(white)) xlabel(30(10)65,labsize(vlarge)) ylabel(0(0.1)0.2,labsize(vlarge))
		graph export "$resultpath/equityshares_housing.eps", replace
		graph export "$resultpath/equityshares_housing.pdf", replace
		
		restore
		
		************** gender gap **************

		// point value
		gen aux = mshare if single == 0 
		bys ageav homeowner: egen aux2 = max(aux)
		
		drop aux
		gen aux = mshare if single == 1
		bys ageav homeowner: egen aux3 = max(aux)
		
		gen diff = (aux2 - aux3)  if homeowner == 1
		replace diff = (aux2 - aux3) if homeowner == 0
		
		// test whether statistically significant
		ttest diff, by(homeowner)
		
		drop aux*
		
		// standard deviation of difference 
		gen aux = std_dev if single == 0
		bys ageav homeowner: egen aux2 = max(aux)
		drop aux 
		
		gen aux = std_dev if single == 1
		bys ageav homeowner: egen aux3 = max(aux)
		drop aux 
		
		gen se_diff = sqrt((aux2)^2 + (aux3)^2)
		
		// confidence bands 
		
		gen diff_high = diff + se_diff*1.96
		gen diff_low  = diff - se_diff*1.96
		
		
		// plotting results
		sort ageav
		
		tw (line diff_high ageav if homeowner == 0, lc(black*0.15) lw(medium)) ///
		   (line diff_low  ageav if homeowner == 0, lc(black*0.15) lw(medium)) ///
		   (line diff_high ageav if homeowner == 1, lc(black*0.15) lw(medium)) ///
		   (line diff_low  ageav if homeowner == 1, lc(black*0.15) lw(medium)) ///
		   (line diff      ageav if homeowner == 0, lp(dash) lc(black) lw(medium)) ///
		   (line diff      ageav if homeowner == 1, lc(black) lw(medium)), ///
		xtitle("Age", siz(vlarge)) ytitle("Gender Gap in %-pts", size(vlarge))  ///
		legend( order(6 "Owners" 5 "Renters" 12 "" 15 "")size(vlarge) position(6) /// 
		col(2) region(lcolor(white)) stack) ///
		graphregion(color(white)) xlabel(30(10)65,labsize(vlarge)) ylabel(-.1(0.2)0.1,labsize(vlarge))
	    graph export "$resultpath/housing_gendergap.eps", replace	
		graph export "$resultpath/housing_gendergap.pdf", replace
		
		drop std_dev diff* se_diff aux*
		
		